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. We study the steady-state motion of mode III cracks propagating on a lattice exhibiting vis- 

' coelastic dynamics. The introduction of a Kelvin viscosity r\ allows for a direct comparison between 

lattice results and continuum treatments. Utilizing both numerical and analytical (Wiener-Hopf) 
techniques, we explore this comparison as a function of the driving displacement A and the number 
^ ' of transverse sites TV. At any TV, the continuum theory misses the lattice-trapping phenomenon; 

/""N , this is well-known, but the introduction of n introduces some new twists. More importantly, for 

large TV even at large A, the standard two-dimensional elastodynamics approach completely misses 
the ^-dependent velocity selection, as this selection disappears completely in the leading order naive 
continuum limit of the lattice problem. 
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Recently, there has been renewed interest in the subject of dynamic fracture [1]. This has been sparked in large 
part by a set of experiments [2,3] that have called into question some of the predictions of the traditional, continuum 
mechanics approach to fracture propagation. Specifically, it has been shown that cracks exhibit a branching instability 
£j ■ long before they reach the predicted limiting speed of advance; this instability leads to enhanced dissipation and 
l — 1 effectively prevents much additional acceleration. Although there are some hints of this instability in the continuum 
approach [4], attempts at a systematic analysis [5] have been inconclusive. 

In this work, we adopt the philosophy of Marder and Gross [6] and consider lattice models of fracture. These 
T^j- i models provide an invaluable test-bed for deciding when and if a continuum formulation is appropriate. After all, 
' if the tip of a brittle crack really occurs at the scale of the lattice, there is no a priori reason for suspecting that a 
7-H , continuum approach could get the correct behavior. It is already clear, for example, that lattice models exhibit a 
sharp (sometimes discontinuous) jump from static cracks to propagating ones; this is not reproducible if one neglects 
lattice scale effects. One might hope, though, that at larger velocity there is some effective continuum description, 
perhaps utilizing the cohesive zone approach of Barenblatt [7]. From our perspective, it is an open issue as to whether 
any such model can accurately predict the behavior of some specific microscopic dynamical system exhibiting fracture. 

Historically, lattice models of fracture received a major impetus from the work of Slepyan [8] , who used the Wiener- 
Hopf technique to solve for steady-state propagation. In his work, he considered the case of infinitesimal dissipation. 
This fact made it difficult to carry out explicit comparisons between lattice results and continuum predictions thereof, 
inasmuch as the latter allows steady-state motion only at the limiting wave speed. One needs dissipative terms so 
as to introduce a new macroscopic velocity scale in order to allow more general steady-state continuum solutions. 
^ ■ Subsequent analyses by Marder and collaborators [6,9] did introduce dissipation in the form of a Stokes term; however, 
they did not explicitly consider the lattice-continuum comparison. In this paper, we introduce dissipation in a different 
form by adding a Kelvin viscosity term to the equation of motion. We will see that the advantage of this choice is 
that the continuum model is in fact an accurate approximation to the lattice dynamics at large enough rjv. Note too 
that from a physics perspective, this type of viscoelasticity appears to have a marked effect on the crack stability [10] 
and is therefore interesting in its own right. 

In this paper, we choose the simplest nonlinear form for the lattice springs, namely that the spring becomes 
completely broken (with no residual force) once it is stretched beyond some threshold. In a future publication, we 
will extend our analysis and results to a more general nonlinear force law. This generalization appears to change very 
little qualitatively with respect to the steady-state problem, although it is crucial in allowing for a direct calculation 
of the linear stability of the propagating fracture. As mentioned above, this stability issue is perhaps the one of 
most immediate relevance as far as the connection to experiments is concerned. But, as we shall see, the steady-state 
problem we solve here offers quite a few subtle and interesting aspects. 

The organization of the paper is as follows. First, we introduce the lattice model, and discuss its basic energetic 
thresholds. We also briefly discuss the static arrested crack solutions. In the next section, we generalize the procedure 
we employed for finding the static solutions to solve for steady-state moving cracks for the case of one row of mass 
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points (N — 1), employing as a foundation the Slepyan ansatz for the form of the discrete steady-state solution. We 
analyze the dependence of the velocity on the driving displacement, studying the effect of the various parameters. 
The most important of these parameters is the viscosity. In the following section, we compare these results to a naive 
continuum limit, finding that whereas this naive continuum limit successfully reproduces the large- velocity regime, it 
fails to account for the nonexistence of solutions below a driving threshold. We then extend our method to arbitrary 
N, again solving for the dependence of the crack velocity on the driving displacement, as a function of the other 
parameters. Subsequently, we again compare our results to those of our naive continuum limit. In particular, we 
focus on the large N limit, demonstrating how in this limit the standard continuum calculation is recovered. We find 
the surprising result that whereas the large-scale structure of the displacement field is almost completely insensitive 
to the crack velocity for velocities less than the wave-speed, the small-scale structure is extremely sensitive. Thus 
the whole velocity selection is completely a function of the lattice-scale dynamics, which no continuum theory can 
reproduce correctly. We conclude with directions for the future and some general observations. 



II. MODEL, ENERGETICS AND STATICS 

We study in this paper the Slepyan lattice model [8] for Mode III fracture, generalized to include Kelvin viscosity. 
The model consists of N infinitely long rows of mass points (with unit mass) coupled horizontally and vertically by 
damped "springs". The bottom and top rows are anchored by "springs" to lines. The system is loaded by extending 
the top row a distance A. The springs all have spring constant 1, except for the bottom row, which has spring constant 
k. All the springs have a viscous damping rj. The bottom springs break when their extension exceeds some threshold 
e. We label the (scalar) displacement of the (i, j) mass from its unstressed equilibrium position as Uij. The equation 
of motion for mj reads 



for j 7^ 1 with u^n+i = A, and 



Ui,i = (1 + 7]-^)(u i+1 j + Ui-i tj + M;, 2 - 3v,ij) - k6(e - + »7^) u *,i • ( 2 ) 

Of particular interest is the case k — 2, which is equivalent to the problem of an up-down symmetric crack, joined at 
the fracture line by springs of strength 1. 

There are a number of important strain thresholds which can be understood from energetics and statics. The first 
is the point at which the uniformly stressed state cracks catastrophically. For our model in which the bottom spring 
has spring constant k and the other vertical springs has spring constant 1, for a (horizontally) uniform state, the 
equilibrium displacements are given by 



< - TT LlA ( 3 ) 
3 Nk + 1 



so that the strain of each spring is kA/ (1 + Nk), except for the bottom 'fc' spring which has strain A/(l + Nk). The 
system will fail catastrophically if the strain of the k spring exceeds e and this gives us our first threshold, 

Au=e(Nk + l) (4) 

The second strain threshold is the Griffith's criterion, which is the A at which the uniformly strained state becomes 
metastable with respect to the cracked state. The energy per column of the cracked state, Uij — A, is just the energy 
to stretch the k spring to cracking, ^e 2 , whereas the energy of the uniformly stressed state is 

kA 2 

(5) 



u 2(Nk + l 

The cracked state is thus energetically favored when A exceeds 



A G = eVNk + 1 (6) 

Note that this is much smaller than Ajj for large N . 

This system is known to possess stationary solutions which represent semi-infinite arrested cracks. For completeness, 
and to begin to build the machinery we will need to treat the moving crack, we briefly outline the solution for this 
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arrested crack. We choose x — to be the position of the last uncracked spring. We solve separately the problem in 
the uncracked (x > 0) and cracked (x < 0) regions and then tie the answers together. To solve, we need to know the 
normal modes of the vertical springs in the two regions. We define the general N x N coupling matrix as 



M N (m) 



-(m + 1) 1 

1 -2 1 

1 -2 1 



1 -2 1 
1 -2 



(7) 



The coupling matrix on the uncracked side is M.^{k) while on the cracked side it is A^at(O). Denote the eigenvectors 
of Mn(0), M.N{k) as £„,E n , with eigenvalue A„, A„, i = 1, . . ., N. (Here and in the following lower (upper) case 
symbols refer to quantities on the cracked (uncracked) side). 
The equation of equilibrium on either side reads 



= U l+1 ,j - 2UiJ + Ui-lJ + MN;3,j'{m)Ui t ji 

The general decaying solution on the uncracked side, i > 0, is 



JV 



= ^ + ^A„(r„) i (s„) J 



(8) 



(9) 



where vF is the uniformly strained solution presented above and 



r„ = 1 - -A„ 



-A„ + i (A„)2 



(10) 



governs the spatial decay of the nth mode, and satisfies (I\) 2 — (2 — Aj)r, + 1 = 0. 
The solution on the cracked side, i < 0, is similar: 
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l ij = A -^2 a n("f n y{€n)j 



(11) 



where 



7r 



= 1--A„ + 



"At. 



r(A„) 2 



(12) 



This solution has 2N unknowns, {A n ,a n }. The equality of the two different expressions for u j provides N 
equations. The equation of motion for x — provides the other TV equations. Solving this 2A^ x 2N inhomogeneous 
system yields the desired answer. The range of validity of this solution is determined by the conditions that tto,i < 
e < it-1,1, so that the spring at i — is the last unbroken spring. 

Doing this, we find that the possible A's span a range, < A < A^ that encompasses the Griffith's value Ac- 
Above A^ the crack has to run, for it has no other alternative. Below A^, any initial crack would head itself. 

In Fig. 1, we look at A^/Ac as a function of width for the natural case, k = 2, and the case k = .2, where 
the material has been weakened along the incipient crack surface. As can be seen, the effect of the width, N, is to 
widen the window, but the effect is quite small, once things are normalized to Ac- The convergence is numerically 
consistent with a 0(1/N) behavior. We see that k, on the other hand, has a dramatic effect, closing the size of the 
allowed window significantly. 

III. MOVING CRACKS, N = 1 

We now look at moving cracks, starting for simplicity with the case N = 1. The key is the Slepyan traveling wave 
ansatz 



Ui(t) = u(t — x/v) 



(13) 
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We thus only need to find the function of one variable u(t). Plugging this ansatz into the equation of motion, we get 
a differential-difference equation which is non-local in time: 

u(t) = (1 + vf f )[u(t + l/v) - 3«(t) + u(t - 1/v)] - k6(e - u(t))(l + j f )u{t) + A (14) 

We choose t — to be the moment at which u exceeds e, so we can replace the step function above by 0(—t). As in 
the static case, we solve the equation separately in the cracked (t > 0) and uncracked (t < 0) regions. It is convenient 
to discretize time with a small time-step dt, so that U — i dt. Then the solution for the uncracked side is 

u(U)=u u +J2M^lY (15) 

I 

where now the T are given by those roots of 

[r n " - 3 - k + T- ,lb ] . (16) 

which lie outside the unit circle. The number rib = l/(vdt) is constrained to be an integer, which implies that our 
resolution in v is limited by our resolution in dt. There are 2nb + 1 roots of this polynomial, (for rj ^ 0), some number, 
n u , of which lie outside the unit circle and thus give rise to a u which converges as t — > — oo. Thus, the solution for 
negative t is parameterized by n u coefficients Ai, I = {1,2, ... ,n u }. Similarly, we solve in the cracked region, and the 
solution is now parameterized by n c coefficients ai corresponding to the roots of Eq. (16) (with k set to zero) which 
lie inside the unit circle. It can also be shown that for sufficiently small dt, n u + n c = 2rib + 1. Thus the entire solution 
is parameterized by 2rib + 1 parameters. As in the static case, the two solutions overlap, this time for 2rib values of t i} 
i = —rib, —rib + 1, . . . , 7i& — 1, so the last uncracked equation at i = —1 fixes u out to i — rib — 1 and the last cracked 
equation at i = 1 likewise fixes u down to i = —rib- The identity of the two expressions for u in the overlap region 
give us 2r%b equations. The last equation we need comes from the equation of motion at the crack point i — 0. Solving 
this inhomogeneous system gives us our desired solution. Reading off ui(0) = e gives us the relation between v and 
A/e we need. 



-, k=0.2 

+. *=0.2 

-k=2 

+, k=2 



I 

1 



1.8 
1.6 
1.4 



v- 2 

1.0 



0.8 - 




100 



FIG. 1. A^/Ag vs. N for the case k = 2, 0.2 



Again, as in the static problem, there is a consistency constraint on the solution, namely that u±(t) not reach e 
before t = 0. In the 77 = problem studied by Marder, this happens for too small v. This holds true, as we shall see, 
for sufficiently small 77. 

In Fig. 2, we present data on v versus A/ Aq for two values of rj = 0.2, 2 with k = 2. Two striking features present 
themselves. First is the divergence of the velocity, which occurs for both values of rj at Ay which we calculated in 
Sec. 2 as the displacement for which the entire system wants to break apart. Note that in this model, there is nothing 
wrong with v > 1, the wave speed in our units. The second important feature of these curves is the different behaviors 
exhibited at the left edge of the graph. The low rj graph exhibits the typical behavior of a subcritical bifurcation 
which ends at a square-root type cusp. The continuation past the cusp to even lower velocities is a numerical artifact 
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of finite dt and vanishes in the dt — > limit. This cusp feature persists in the 77 — > limit, and for smaller velocities 
the solutions are inconsistent with the condition mentioned above that u(i) < e for t < and so are unphysical. For 
larger -q, where the system is sufficiently overdamped, the solutions persist to zero velocity. However, the dependence 
on A/Ag is very singular, and the graph approaches zero velocity at A^ exponentially in l/v. Before proceeding to 
a full survey of parameter space, it is useful to develop a naive continuum limit for our system, which is analytically 
more tractable and serves as a useful benchmark for our discussions. We do this in the next section. 
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FIG. 2. v versus A/ Ac for r\ — 0.2, 2 for N = 1, k = 2. The calculation was done with dt — .05. The solid lines are the 
naive continuum results for the same parameters. 

IV. NAIVE CONTINUUM LIMIT 

We now develop the naive continuum limit of the equation of motion for N = 1 steady-state moving cracks. We 
obtain this limit by simply replacing the finite difference by a derivative, yielding: 

The solution for t < is 

u(t) = I ^+A 1 e^ t (18) 

where Q\ is the unique root with positive real part of the polynomial P (— (1 + k); Qi) where in general P(\; Q) is 
defined by 

P(A; Q) = V vQ 3 + (1 - v 2 )Q 2 + (1 + V vQ)X . (19) 

In general, for A < 0, P has one positive root, with the other roots having negative real parts. Similarly, the solution 
for t > is 

u{t) = A - a 2 e- q2Vt - a 3 e- q3Vt (20) 

where — <j2, — <Z3 are the two roots with negative real part of P(— 1; q) One can solve for A\, a 2l and 03 by using the 
continuity of u and its first two derivatives at t = 0. 

One can learn a few things from these equations. As pointed out by Langcr in his study of a related model [11], 
vanishing 77 in the continuum is a singular limit, since it controls the highest derivative. Secondly, the special role of 
the wave speed, v = 1, is clear. Most interesting, however, is the question of when we expect this continuum limit 
to be valid. The condition is that at least some of the exponential decay rates are small, such that the solution will 
look smooth on the lattice scale. For small v, Qi <~ \/k + 1, q 2 ~ 1, and q% <~ l/(i]v), so none of the g's are small 
and the continuum limit is not reliable. For large v, things are different. Here, Qi <~ v/rj, g 2 ,3 ~ (f] ± y/rf — &)/(2v). 
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The large value of Q\ corresponds to the existence of a boundary layer which allows for the matching of the highest 
derivative term. If this boundary layer does not affect the lower order matches (as is fairly typical), then we would 
expect that the continuum result will agree with the lattice answer. 

It is straightforward to work out the small velocity limit of the continuum theory. Using the limiting values of 
the q's above, and solving the linear system, one gets that e <~ A/\/fc + 1 + 0(v), so that A <~ Ac + 0(v). Thus, 
the continuum solution starts at the Griffith's point Aq with zero velocity and the velocity grows linearly as A is 
increased. Continuing the calculation to first order in v, we find v = (A/Aq — l)/[(y/k + 1 — 1)77]. Thus, the velocity 
is inversely proportional to 77. This is as expected, since for v -C 1, the velocity only enters in the combination r/v. 

The large velocity limit is also analyzable. Again using the g's found above, we obtain to leading order that 

v = [(1 — A;y)/(fcr7 2 )] ^ 4 so that v diverges at Au- Note also that in this regime v scales as y/rj, so that as viscosity 
increases so does the velocity. This is because the only reason that propagation at v > 1 is possible is because of the 
viscosity, so the larger the viscosity the more efficient the propagation. 

As 77 goes to zero, the continuum limit must break down. The velocity increases very rapidly (at a rate proportional 
to I/77) to near 1, and stays there till A is near Au, whence it rapidly diverges. The velocity crosses unity at 
A = (1 + fc) 2 / 3 with slope of order rj 2 ^ . Thus, with vanishingly small 77, steady-state propagation, at least in our 
naive continuum limit, is only possible at the wave speed v = 1. 

We are now in a position to compare our lattice calculation to the continuum results. For example, Fig. 2 above 
also displays the continuum curves. We see that in general the continuum calculation is very good only at the largest 
velocities. The agreement gets progressively worse for smaller velocities and breaks down completely for A less than 
the arrest value A^. This must be since the continuum calculation has the velocity going smoothly to zero at Aq, 
which it never does due to the existence of arrested solutions. Also, the agreement is better at larger 77; above some 
critical 77, (about 77 m 0.5 for k = 2), the continuum curve serves as an upper bound on the lattice curve and thus 
is a fair approximation until we start getting close to the arrest value. All of this is to be contrasted to what would 
have been obtained had we set 77 = and introduced instead a Stokes velocity term proportional to u. Then, the 
continuum theory has a limiting velocity of v = 1 and since this feature is not shared by the lattice dynamics, the 
continuum approximation would be nowhere accurate. 
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FIG. 3. r]v versus A/Ag for 77 = 1, 1.5, 2, 4, 8 for N = 1, k = 2. The calculation was done with dt = .05. 

We saw from the continuum calculation that for velocities v <C 1, the relevant parameter for the continuum 
calculation is 777;. We can test this for the lattice model by plotting 777; versus A. This is presented in Fig. 3. We see 
that asymptotically for large 77 this scaling sets in, but for finite 77, the existence of the lattice length scale ruins the 
simple scaling. 

It is interesting to see what happens for smaller k. We saw that the window of arrested solutions is significantly 
smaller for smaller k. In Fig. 4, we present the analog of Fig. 2, but this time with k = 0.2. We see that the results 
are as far from the continuum limit as they were with the larger fc. While the window of arrested solutions is smaller, 
so is the value of Au/ 'Aq, so the entire picture is just shrunk to a smaller range of A. 
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A. Wiener-Hopf Solution 



Before we turn to general N, we will present the equivalent Wiener-Hopf solution of the N = 1 problem. This 
method is more involved for the case N = 1, but it is a model for the Wiener-Hopf solution we will present for general 
N, which we allow us to draw analytical conclusions in the large N limit. The basic method follows that of Marder 
and Gross [6], but the presence of viscosity adds some new twists which are worthy of comment. 
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FIG. 4. v versus A/ Ac for r\ — 0.2, 2 for N = 1, k = 0.2. The calculation was done with dt = .05. The solid lines are the 
naive continuum results for the same parameters. 

To begin, we define the Fourier-transform of the right- and left-hand pieces of the u field as follows: 

/oo 
vdt 9(±t)e iKvt u(t) (21) 
-OO 



u ± ' 



It should be noted that are analytic in the upper and lower half-planes respectively. The Fourier transform of the 
u field, u, is just the sum of the two parts: u{K) = u + + vT . In terms of these fields, the equation of motion reads 



= (iijvK 3 - (1 - v 2 )K 2 ) u - (1 - irjvK)u - fc(l - ir)vK)u~ + AS{K) - krjvu{0) 



(22) 



The last term is noteworthy, and arises because the time derivative does not act on the ^-function in the last term 
in Eq. 17. Expressing u in terms of its component pieces, we recognize that the coefficients of are nothing put 
P(— 1; —%K), P{— (1 + fc); —iK) respectively, that we encountered in our solution above, so that the the equation of 
motion reads 



= P(-(l + k); -iK)u~ + P(-l; -iK)u+ + A5(K) - kr]vu(0) 

Next, we factor the P's in terms of their roots 

P(-(l + fc); -iK) = ir)v(K - iQ x ){K + iQ 2 )(K + iQ 3 ) 
P(-l; -iK) = irjv{K - iq x ){K + iq 2 )(K + iq 3 ) 



Dividing through by ir\v(K — iq±)(K + iQ2){K + 1Q3) we obtain 

A 



K-iQ x __ (K + iq 2 ){K + iq 3 ) 
U = —r, : U + T~rz — ; . - ; U 



K - iqi 



(K + iQ 2 )(K + iQ 3 ) 



rjvqiQzQi 



S(K) 



ik 



(K-i qi )(K + iQ 2 )(K + iQi 



MO) 



where we have used the S(K) to simplify its prefactor. We use the fact that q\q 2 qz = — P(— 1; 0)/(r)v) 
rewrite our equation as 



= 



K-iQx„_ 

1? : — u 

K - iq 1 



(K + iq 2 )(K + iq 3 ) „ + Aq 2 q 3 

-u — - - o(K) 



ik 



(K + iQ 2 ){K + iQ 3 ) 



Q2Q3 



(K -iqi)(K + iQ 2 )(K + iQ 3 ) 



*(0) 



(23) 

(24) 

(25) 
1/rjv to 

(26) 
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To proceed, we have to decompose the two inhomogeneous terms into pieces analytic in either the upper- or lower- 
half-planes. The result is 

n K-iQ x __ (K + iq 2 ){K + iq 3 ) „ + Aq 2 q 3 
= — : U + —— — , . T , tt^-U — 



K-i qi {K + iQ 2 )(K + iQ 3 ) Q 2 Q 3 \K + iO+ K — i0+ 

k ( 1 K + i{ qi + Q 2 + Q 3 ) 



(qi + Q 2 ){qi + Qs) \ (K- %) (K + iQ 2 ){K + iQ 3 ) 



u(0) (27) 



The key to the Wiener-Hopf method is the realization that the sum of the terms analytic in either half-plane have to 
vanish, allowing us to solve for u ± . We find 

.+ = iAg 2 g 3 (K + iQ 2 )(K + iQ 3 ) ik(K + i( qi +Q 2 + Q 3 )) 

U ~ Q2Qs(K + iO+){K + iq 2 ){K + iq 3 ) ( Ql + Q 2 )( Ql + Q 3 )(K + iq 2 )(K + iq 3 ) H ' 

— = -iAQi(K - i qi ) ik 

(l + k)q 1 (K-iO+)(K-iQ 1 ) + (q 1 +Q 2 )(q 1 +Q 3 )(K^iQ 1 ) U[ > { ' 

Notice that the poles in give rise to exactly the same exponential terms in u that we found previously. It can be 
explicitly verified that the two forms of the solution are equivalent. For our purposes, it is sufficient to examine what 
happens in the small v limit. Then, as we have already noted q 3l Q 3 ~ l/{riv). The first, A, term on the right-hand 
sides approaches a finite limit, with 0(v 2 ) corrections, whereas the second ^(0) term vanishes linearly in i]v. More 
explicitly, we find 

••• - i^(K + iQ 2 ) //."/„■ _ u(0) (29) 



QiQsiK + iO+)(K + iq 2 ) ( qi + Q 2 )(K + iq 2 ) 



Using q\ = q 2 = 1, Qi = Qi = a/1 + k, and q 3 = Q 3 = l/rjv, and using the inverse Fourier transform to evaluate this 
in the limit x — > + , we obtain 

u(0) - -JL= kvV u(Q) (30) 



Recognizing u(0) = e and vT+fc = Ag/e, after reorganizing we obtain 



/A \ i + VT+k 
^ V = \ Ag _ V k (31) 

which is easily seen to be equivalent to what we obtained from the direct method. As can be appreciated, the Wiener- 
Hopf method is much more involved than the direct method. Nevertheless, it will be the essential tool for analyzing 
the large- N limit. 



V. GENERAL N 



It is straightforward to extend the lattice calculations to arbitrary N. The basic method is the same: we solve the 
problem on the two sides of the crack tip position and patch the two solutions together. The solution on either side 
is again a sum over modes, which are a direct product of modes in the vertical direction, given by the eigenmodes of 
A^Ar(m = 0, k), with the modes in the horizontal direction. Thus there are a total of Nn u and Nn c modes on the 
uncracked and cracked sides, respectively. The solutions on either side have to overlap for each value of the vertical 
component j, so there are an appropriate number of equations for the unknown coefficients of each mode. As for 
N = 1, the condition u(0, 1) = e is used to determine the driving A corresponding to a given velocity. 

We can also generalize our continuum calculation to finite N. As in the N = 1 case, we replace finite differences in 
t with derivatives, giving us N coupled third-order differential equations. Again, we can solve these exactly on either 
side of the crack tip t — 0, and match the functions and their first and second derivatives at this point. The functions 
are characterized by N modes on the uncracked side, with decay rates <3i, n , an d modes on the cracked side, with 
decay rates q 2 , n , 93,n- For each n, Q\, n is the positive root of of the polynomial P(A„) defined in Eq. (19) above. Let 
us denote the other roots of this polynomial, which we will need later, by —Q 2 _ n , —Qs t n- Similarly, —q 2t n, — q 3 _ n are 
the two negative (real part) roots of P(A„). with the third, positive, root being labeled by qi tTl . Implementing these 
procedures, we again calculate the crack velocity as a function of A/Aq. Again, we compare this data to that of our 
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naive continuum (in x) calculation for the same value of N. We present in Fig. 5 the results for our overdampcd 
case rj = 2, for N = 1,2,5,10,15. Qualitatively, not much changes with N. The most important feature in that 
the middle section of the curves get progressively flatter as N increases. This must be the case, since the point of 
divergence, Ay, measured in terms of Aq, increases as iV 1 / 2 . We also note that the data for low velocities seems to 
converge fairly rapidly as N increases, and the rate of convergence slows as v increases. Again, as in the N = 1 case, 
the continuum results accurately reproduce the lattice calculations for large v and are completely wrong for small v, 
missing the lattice-induced arrest phenomenon. 
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FIG. 5. v versus A/ Ac for r\ = 2 for N = 1, 2, 5, 10, 15; k — 2. The calculation was done with dt — .1. The solid lines are 
the naive continuum results for the same parameters. 



VI. LARGE N 



The physical problem of cracking a macroscopic object corresponds to the limit of large, but finite, N. The lattice 
calculations are prohibitively expensive for too large N . However, our naive continuum calculations can be carried 
out for fairly large iV's. Using the fact that for small v, the convergence in N is rapid, and for larger v, the naive 
continuum results are reliable, we can piece together a fairly complete picture of what we expect in the macroscopic 
limit. In particular, it is interesting to compare this with the standard continuum theory in order to understand the 
limitations and successes of that theory. 
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FIG. 6. v versus A/Ag in the continuum approximation for for N = 10, 50, 200, 400, with r) = 2, k = 2. 
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To begin, we present in Fig. 6 the results of our naive continuum theory, extended to larger values of N. The most 
striking feature of this graph is the slow convergence that sets in near v — 1. Exploring numerically, we find that for 
fixed v < 1, the data converges with large N as N~ x l 2 . However, the coefficient of this N~ x / 2 correction becomes 
ever larger as v approaches unity. Looking at the value of A/ Ac where v = 1, it appears to be diverging as TV 1 / 6 as 
N — > oo. Thus, in the macroscopic limit, the crack speed is effectively bounded by the wave speed. 

To proceed further in studying our naive continuum theory at large N, it is useful to derive the Weiner-Hopf 
solution. To do this, we first Fourier transform the fields, writing 

/°° dk 
^_ e -ivKt~. {K) (32) 

The equations for all j ^ 1 are translationally invariant in t, and become algebraic. The structure of these equations 
is 

= (ir/vK 3 - (1 - v 2 )K 2 ) Uj + (1 - ir]Kv)M N I \Uj> + 5 ji2 (l - ir]vK)u x + 5 hN AS(K) (33) 

Defining f{K) = (irjvK 3 — (1 — v 2 )K 2 )/(l — irjvK) and denoting the n x n identity matrix as l n , we can using 
Cramer's rule explicitly solve for u 2 in terms of u\ as follows 

_ det N _ 2 (f(K)l + M(l)) (-l) jV A 
U2 ~ -det N ^(f(K)I + M(l)) Ul ~ det N ^M(l) 6{K) ' (34) 

where in the last term we have used the 5(K) to simplify the determinant. To treat the j = 1 with its step functions 
we define 

/oo 
vdt 9(±ty Kvt Ul (t) (35) 
-oo 

so that u\ = u + + u~ . The j = 1 equation now reads 



= (1 - irjvK) 



tf lK\ 1V det N _ 2 (f(K)l + M(l)) „ 



(36) 



-W°)- d ^MW S{K) (3?) 
Multiplying through by det N - X (f(K)l + M(l)) /(l - i-qvK), we get 

= det N (f(K)l + M{0))u- + det N (f(K)l + M{k))u+ 

l^iryvK M{K) (38) 

The determinants are easy to calculate in the diagonal bases of the M's, and have zeros at K's corresponding precisely 
to i times the roots of the polynomials P(A„), P(X n ) we encountered in our original real-space calculation. We can 
thus write 

N 

-N 



det N (f{K)l + M(k)) = (1 - i V vK))- N 11 P{K n --iK) 

71=1 

, . vJV JV 

= [ YZJ^k ) Yl( K -iQi,n){K + iQ 2>n )(K + iQ 3tN ) (39) 



n=l 

and similarly 

N 



det N (f(K)X + M(0)) = (1 - ir]vK)- N J] P(A„; -%K) 

n=l 

, . -.NN 

= L^vK ) Il( K - i ^)(K + iq 2 , n )(K + iq 3)N ) (40) 

^ ' n=l 
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Similarly, if we denote the eigenvalues of .Mjv-i(l) by £ m , m = 1, . . . , N — 1, we can express detN-i(.f(K)T + M(l)) 
in terms of the roots xi,m, ~X2,m, and -X3,m of P(^ m ) 



AT-1 



de^_ 1 (/(X)X + M(l)) = (1 - i^vK)- 1 *- 1 J] P(£ m] -iK) 

irjv \ 



JV-l iV-l 



m— 1 



We can then re-express Eq. 38 as 

^TlmC^ ~ ^Xl,m)(^ + iX2.m)(K + %,m) ,^ A 

n„(^-%,n)(^ + ^2,„)(# + iQ 3 ,n) " U ' M W n„9l,n<32,„Q3,n 1 ' 

n K-iQ^n — TT (-K" + ig2,n)(^ + Zg 3 ,n) -+ f,^ 
a ^ (i^ + ig 2 ,n)(^ + zQ3,n) 1 J 

Un( K -m,n)(K + iQ2,n)(K + iQ 3 , n ) lW 11 Q 2 ,„Q 3 ,n 1 ' 1 j 

where in the last step, we applied the identity 

I] «l, n «2,„«3,n = (-w)- N detM N {0) - • (45) 

n 

and where we have used the easily verified fact that detJVl^(m) = (—l) N (mN + 1). Note that this product result 
nicely reduces to the result we previously obtained for N = 1. Again, as in the N = 1 case, to proceed with the 
Wiencr-Hopf method, we need to break up the last two terms into pieces analytic in the upper and lower half-planes. 
The ui(0) piece does not appear to have a simple breakup. However, for large N, the effect of this term becomes 
irrelevant, since wi(0) is a factor N 1 / 2 smaller than A. The last term is easily broken up as we did in the N = 1 case. 
We find that to leading order 

(K + iQ 2 , n )(K + iQ 3in ) 
'\K + iO+J \i-Q 2tn Q^ n (K + iq 2in )(K + iq^ n ) ■ 1 > 

In the large TV limit, we can use this to evaluate u explicitly. (We need not concern ourselves with u~ , since for 
t < 0, u\(t) is always smaller that ui(0), and so does not contribute to leading order.) To proceed, we need the explicit 
form of the q, Q's to leading order. As we shall see, the behavior is controlled by modes where n « N. For these 
modes, we may approximate the eigenvalues of .Mjv(fc) by Ai ; „ = A 2 , n = —n 2 ir 2 /N 2 7 so that Qi^ n = Q 2 ,n = N ^/i_ v -i > 
Q 3 ,n = (1 — v 2 )/(i]v). (Here we have assumed that v is less than and not too close to 1.) Similarly, for A^at(O), wc find 
Ai,„ = \ 2 , n = -(n - ±) 2 tt 2 /N 2 , so that q 1<n = q 2 , n = j^=r, q 3 ,n = (1 - v 2 )/(tjv). Notice that since Q 3>n w q 3>n , 
the factors involving these quantities cancel. This has the immediate consequence that the viscosity 77 has completely 
dropped out of the problem in this limit. 

The remaining expression has poles at — i0 + and at —iq 2 . n . We can evaluate the residue of each of these poles 
explicitly. The residue at — i0 + is immediately seen to be 

Res(e- lKvt u+)\_ m = 1A . (47) 

Evaluating the residues at the other poles is more complicated. To proceed, let us cut-off the product at some large 
n = N c « N. Then, our approximate expressions for the q's and Q's are valid. This leads to 

iMe-™ 4 +) Ul ,„ . w + w* - iw. + 1 - ») 



«+ 



(n 



(n - ^)T(N C + l)r3(I)r(n)r(A c + 1 - n) 



-** N ^ . t\ i' ( \ + °(jr ) ( 48 ) 
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We then take the limit of N c — > oo to find the final answer for the macroscopic displacement field 




12 3 

vt/N 

FIG. 7. wi(t)/A versus vt/N in the continuum approximation for k = 2, 77 = 2, TV = 25, compared with the large- TV analytic 
result, Eq. (49). 

This final answer exhibits the well-known square root branch cut at the crack tip location, t = 0. It is worth noting 
that this behavior of the displacement gives rise to a macroscopic stress field which actually diverges as ^375 (recall 
the extra derivative due to the Kelvin viscosity) near the crack tip. This surprising finding renders invalid the 2-d 
continuum calculation of Langer [11] who studied this problem with the additional complication of a finite length 
cohesive zone. A correct continuum formulation which does reproduce the essential formula Eq. (46) is presented in 
the Appendix. 
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FIG. 8. u\(t) / UAsymptit) versus vt in the continuum approximation where UAsympt is the large- TV analytic result, Eq. (49), 
for TV = 1, 2, 5, 10, 25. Again, k = 2, 77 = 2. 

A comparison of the above prediction with the numerically computed displacement is shown in Fig. (7) for the 
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case N = 25, plotted as a function of the macroscopic scale vt/N. We see that our large N analytic result correctly 
reproduces the large-scale structure of the crack displacement. It does worse, though still quite accurate, close in to 
the crack tip. To demonstrate this more quantitatively, we present in Fig. (8) the ratio of the crack displacement 
u\(t) to the large- N analytic result, for various N. Now the data is plotted as a function of the microscopic scale 
vt. We see that curves are all quite similar. They have a square-root divergence at the origin, since the analytic 
prediction is that u(t) vanishes at t = 0, whereas the true answer is finite. By N = 25, they have converged to a 
limiting curve. This means that the finite- N theory possesses a well-defined microscopic structure, in addition to the 
universal macroscopic structure defined by the standard continuum theory. This microscopic structure is on the scale 
of the lattice constant (in the y-dircction) and is of course invisible to the standard continuum theory. 

This observation implies that the entire issue of velocity selection via the condition that u(x = 0, y = 1) be fixed 
to equal the breaking displacement is out of reach of the leading order macroscopic limit. Thus, as an example, 
the velocity depends explicitly on r\ (as opposed to the macroscopic displacement which explicitly does not) even 
for arbitrarily large N. Conversely, calculating the macroscopic field in the continuum limit does not suffice for 
determining the crack speed which is always fixed at the lattice scale. Situations with equivalent macroscopic fields 
can have arbitrarily different crack velocities. 

VII. SUMMARY 

In this paper, we have studied in some detail the steady-state motion of mode III viscoelastic cracks in a lattice 
model of the microscopic dynamics. The most important finding are 

1. The existence of a minimum velocity for crack propagation is dependent on the viscosity. At low rj (and indeed 
in the lattice models without dissipation that have been studied to date), the steady-state branch starts at finite 
v. For highly damped systems, on the other hand, the branch extends all the way to v — at the upper end of 
the allowed range of A for arrested cracks. 

2. For finite N, a continuum approach (in x) for the crack does accurately predict the lattice results for values of 
the driving away from the lattice trapping (low or zero velocity) regime. 

3. Taking the macroscopic limit (N — ► oo) allows us to recover the expected macroscopic behavior that the 
displacement grows as ^fx — a^ip once we leave an inner core region of the lattice scale. The coefficient of this 
term can be calculated by using a continuum theory with the proper boundary conditions. A key feature of this 
macroscopic theory is that the viscosity becomes irrelevant. 

4. However, the velocity selection as a function of the imposed displacement is wholly controlled by the core and 
cannot be accurately arrived at by any theory which does not explicitly consider the lattice scale. In particular, 
viscosity plays a crucial role in this feature of the physics. 

As mentioned in the Introduction, the next step in our research program will be to consider the modifications 
introduced into the aforementioned results by having a continuous but nonlinear force law. In particular, having a 
force which immediately drops to zero means that there is no way that the system could dynamically decide to create 
a "cohesive" zone of mesoscopic (i.e., scaling as a positive power of N) proportions. In such a zone the displacement 
would be such that the force law would be beyond the linear spring regime but not so large that the force would be 
effectively zero). If it were of large size, it would lead to a more complex continuum theory along the lines suggested 
by Langer and co-workers [11,12]; if it were purely on the lattice scale, it would change nothing. Our initial evidence 
suggests the latter, and we hope to report on this in the near future. 

As far as the physics of fracture is concerned, we must address several issues that go well beyond the studies in this 
paper. Since most of the experiments concern mode I cracks, wc need to extend our results to that situation; this 
is technically challenging but should not lead to any significant surprises. Next, we must explicitly investigate the 
stability of our steady-state equations. Finally, all lattice models leave out the possibility of ductile behavior involving 
the emission of dislocations from the crack tip; comparison to experiment and to molecular dynamics simulations will 
enable us to learn when these additional phenomena are crucial or alternatively when one can get by with a purely 
"brittle" model. 
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APPENDIX: THE DIRECT CONTINUUM CALCULATION 



In this appendix, we present a direct continuum (in x and y) calculation of the steady-state crack. We will see that 
it recovers directly the leading-order results of the large- N limit calculation presented in section VI above. 
To begin, we write the displacement field u(x, y, t) — u(t — x/v, y) in Fourier space: 

u(t,y) ,Al+f ^ HK ^K Vt sinh(k(W -y) 

where k y satisfies the dispersion relation 

{l-ir,vK){-K 2 + k 2 y )+v 2 K 2 = (A2) 

The crack is chosen to begin at x = so u(x < 0, 0) = 0. On the crack surface y — 0, x > 0, we must set du/dy = 0. 
Note that this condition implies that the normal stress on the free surface, (1 + vanishes. However, it is 

incorrect to assume, as Langer [11] did in a parallel calculation, that the vanishing of the normal stress is a sufficient 
condition, as this allows for (unphysical) displacement fields that do not have du/dy = 0. As we have seen, the 
macroscopic field possesses a square-root singularity at the origin, while Langcr's condition eventually results in a 
much weaker x 3 ^ 2 singularity (in the absence of Barenblatt type surface stresses). Our condition implies 



J — ( 



^u(K)(-k y ) coth(k v W) = ~ + 0{-t)G{t) (A3) 



or, Fourier-transforming this equation: 

u{K){-kyW) coih{k y W) = ^K K ) + G~{K) (A4) 
where G~ is the transform of 0(—t)G and has no zeros or roots in the lower-half-plane. To proceed, we use the identity 

k y W coth{k y W) = n ^f^li (A5) 

Now, we can use the dispersion relation to eliminate k y in favor of K. If we define A„ = — [(n — ^)tt/W] 2 and 
A„ = — (nn/W) 2 , then we find that 

n— 1 v ' 

Notice that A„, A„ are precisely the same as those in the finite- N calculation for n <C N, if we identify W = N. 
Expressing the P's in terms of their roots, we get 

kyWCOth(kyW) - qinq2nq3n{K _ iQin){K + iQ2n){R + ^ (A7) 

Plugging this into Eq. (A4), and reorganizing, we obtain 

n Q2,nQ3.n( K + i( l2,n)(K + iq 3 . n ) . _ -|-|- qi. n (K - lQi. n ) = ^_ -pr q^njK^-iQj^n) 

n q2, n qsAK + iQ2, n )(K + iQ 3 , n ) *■} Qi, n (K - iqi,n) [ ' *■} Qi,n(K - iqi,n) [ ' 

Decomposing the 5-function as in the finite- A'' case, and separating out the pieces analytic in the upper-half plane, 
we get 
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i+ -Q Q2,nQzA K + i( l2,n){K + ig 3 ,n) _ A_ . 



Q2,nQ3,n 

so that 



(K + iQ 2 , n ){K + iQ 3 , n ) 



n ^2,ra</3,n^ t l^njjn T <-W3,n) (AW) 
a Q2,nQ3,n( K + il2,n)(K + iq 3 , n ) 1 ' 

That this result is the direct equivalent of our leading-order finite- N result, Eq. (46), is clear. One word of in- 
terpretation is called for, however. To achieve the macroscopic limit of our finite-iV result, we needed to take the 
width N large. This in turn implied that viscosity was irrelevant in the macroscopic limit (unless we scaled it by 
a power of N with no obvious physics justification). If we work directly in the continuum, however, we obtain the 
same final result without having to take W large. Thus, the ratio of the viscous length scale to W is arbitrary in 
this continuum calculation. Nevertheless, if we examine the large- W limit of our continuum calculation, we again will 
find that viscosity becomes irrelevant. It is also worth reiterating that this continuum calculation has no sign of the 
subdominant pieces which control velocity selection. 
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